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We perform molecular dynamics simulations to study lithium dynamics in a model of LiP03 glass 
at temperatures below the glass transition. A straightforward analysis of the ionic trajectories shows 
that lithium diffusion results from jumps between sites that are basically unmodified on the time 
scale of the lithium ionic relaxation. This allows us a detailed identification and characterization 
of the sites. The results indicate that the number of lithium sites is only slightly bigger than the 
number of lithium ions so that the fraction of vacant sites is very limited at every instant. Mapping 
the ionic trajectories onto series of jumps between the sites provides direct access to lithium jump 
dynamics. For each site, we determine the mean residence time r s of an ion and the probability 
Pa that a jump from this site to another site is followed by a direct backjump. While a broad 
distribution G(lgT 8 ) shows that different sites feature very diverse lithium dynamics, high values of 
pi give direct evidence for correlated back-and-forth jumps. A strong decrease of p h a with increasing 
r a indicates that the backjump probability depends on the dynamical state of an ion. Specifically, we 
find that correlated back-and-forth jumps are important at short times in the relaxation process, but 
not on the time scale of the lithium relaxation, where the hopping motion resembles a random walk. 
We further study how the local glass structure and the local energy landscape affect lithium jump 
dynamics. We observe substantial effects due to the energy landscape, which are difficult to capture 
within single-particle approaches. Our results rather imply that lithium migration is governed by 
the competition of the ions for a small fraction of vacant sites in a disordered energy landscape. 
Consistently, a statistical analysis shows that a vacancy mechanism dominates the repopulation of 
the lithium sites. 

PACS numbers: 66.30.Dn 



I. INTRODUCTION 

The mechanism for fast ion transport in glasses has 
challenged the combined efforts of industrial and scien- 
tific research over the last decades It is now well es- 
tablished that the macroscopic charge transport results 
from hopping of mobile ions between localized states in 
the glassy matrix, 1,3 However, the mechanism of this hop- 
ping motion is complex and several coexisting models 
seek to rationalize the very rich phenomenology of ionic 
glasses. While several authors postulated that ionic mi- 
gration is governed by the Coulomb interactions among 
the mobile ions^ 1 ^ others stressed the relevance of the 
disordered energy landscape provided by the glassy ma- 
trix^SiS and percolation approaches were applied for 
a quantitative description of ionic diffusion>i2*ii*i2i The 
combined effect of these interactions was also regarded 
as important In addition, it was argued that 

structural properties, e.g., microsegregation of the mo- 
bile ions, play a major role for ion dynamics i 16 : 17 Also, 
the relevance of substantial relaxation of the glass struc- 
ture below the glass transition temperature T g was em- 
phasized^ 

To decide between different models of ion transport in 
glasses, detailed information about the ionic motion is 
required. For example, it is interesting whether all ions 
make the same contribution to the macroscopic charge 
transport or whether a distribution of mobilities exists. 
Further, the relevance of correlated back-and-forth jumps 



needs to be quantified. Due to this effect, which is one 
of the cornerstones of various modelling approaches, the 
dc conductivity Ode can be much smaller than expected 
based on the ionic jump rates. These points are related 
to the question about the origin of the nonexponential 
ionic relaxation, observed in electrical and mechanical 
relaxation studiesi2ii2i22i2ii22 Such behavior can result 
from two fundamentally different scenarios S In the het- 
erogeneous scenario, all particles are random walkers, 
but a distribution of jump rates exists. In the homo- 
geneous scenario, all particles obey the same relaxation 
function that is intrinsically nonexponential due to corre- 
lated back-and-forth jumps. Additionally, it is intriguing 
to study how ion dynamics depends on the local glass 
structure. 

Measuring the frequency dependent electric conduc- 
tivity <j(v), an experimental study of ion dynamics is 
possible. While the frequency independent value at low 
frequencies, a da is an important material parameter, 
a strong dispersion at higher frequencies indicates the 
presence of correlated back-and-forth motions ASiS^SiSi 
Multidimensional nuclear magnetic resonance (NMR) ex- 
periments provide direct access to ionic jump dynam- 
ic i 28 4 29 J 30 4 3i i 32 Applications on crystalline and glassy sil- 
ver ion conductors showed that the nonexponential de- 
population of the silver sites is due to a distribution of 
jump rates, i.e., dynamic heterogeneity, rather than to 
an intrinsic noneponentiality^i^S In other words, no ev- 
idence for back-and-forth jumps was found, at variance 
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with conclusions from electric conductivity studies. 

In molecular dynamics (MD) simulations, the trajec- 
tories of all particles are known and, thus, unique mi- 
croscopic insights are available. Applications on ionic 
glasses confirmed that the dynamics of the mobile ions 
can be decomposed into vibrations about sites and occa- 
sional jumps between them, 33 i 34 i 35 i 36 i 37 i 38 i 39 i 40 i 41 i 42 i 43 i 44 
Most simulation studies of ionic jump dynamics focused 
on alkali silicate glasses. For lithium silicate glasses, it 
was demonstrated that both homogeneous and hetero- 
geneous dynamics contribute to the nonexponentiality 
of ionic relaxationi 4 ^ 4 ^ For sodium silicate glasses, no 
evidence for back-and-forth jumps was found, but there 
are dynamic heterogeneitiesi^i 4 ^*^ 4 ^^ Specifically, fast 
ion transport was observed along preferential pathways, 
where, however, the presence of these channels was not 
attributed to a microsegregation of the sodium ions, as 
proposed by other workers^&^S 

Very recently, we performed MD simulations to study 
ion dynamics in an alkali phosphate glass44 Based on 
multitime correlation functions, we demonstrated that 
the nonexponential lithium relaxation in LiPC>3 glass re- 
sults from both correlated back-and-forth jumps and a 
broad distribution of jump rates. A quantitative anal- 
ysis showed that the relevance of the heterogenous con- 
tribution increases with decreasing T. In addition, we 
observed an exchange between fast and slow lithium ions 
that occurs on the timescale of the jumps themselves. 
Thus, the dynamic heterogeneities are short lived, indi- 
cating that sites featuring fast and slow lithium dynam- 
ics, respectively, are intimately mixed. 

Heuer and coworkers^ used a computational approach 
to demonstrate that determination of the lithium sites 
in lithium silicate glasses provides new insights into the 
mechanism for ionic diffusion. A complete identification 
of the sites from the ionic trajectories showed that the 
number of sites is only slightly bigger than the number of 
ions, implying that ion dynamics is most appropriately 
described in terms of mobile vacancies. For a detailed 
study of the jump diffusion mechanism, the alkali tra- 
jectories were mapped onto sequences of jumps between 
the sites^i An analysis of these sequences revealed that 
the tendency for back-and-forth jumps depends on the 
dynamical state of an ion. Specifically, the back-jump 
probabilities are high when the lithium ions have escaped 
from sites characterized by comparatively short residence 
times. 

In the present MD simulation approach, we investigate 
lithium jump dynamics in LiP03 glass via identification 
of the lithium sites. In particular, we study the effects of 
the local glass structure and the local energy landscape 
on the lithium hopping motion. To unravel the origin 
of the pronounced dynamic heterogeneities observed for 
this model glass in our previous work^ we analyze which 
factors determine the residence time at a site. Such anal- 
ysis is also performed for the correlated back-and-forth 
jumps so that the present study may shed a light on the 
question why previous experimental and computational 



studies came to different conclusions about the relevance 
of this phenomenon. Furthermore, we use a statistical 
approach to ascertain the dominant mechanism for the 
rcpopulation of the lithium sites. All analyses are car- 
ried out in a temperature range, where the phosphate 
glass matrix, apart from local fluctuations, is rigid on 
the 20 ns-time scale of our simulation, while the lithium 
ionic subsystem can still be equilibrated. 

II. MODEL AND SIMULATION DETAILS 

The interactions of the ions in the studied model of 
LiPC>3 glass are described by the potential 
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S o/3 (r)=^^+4tf exp(-r/p), (1) 

where e is the elementary charge and r denotes the dis- 
tance between two ions of type a and (3, respectively 
(a, (3 € {Li, 0,P}). Karthikeyan et alS adjusted the 
parameters of this potential to enable a realistic descrip- 
tion of the structure of LiPC>3 glass. They showed that, 
though there may be some deviations in the intermedi- 
ate range order^S the simulated glass consists of well 
defined phosphate tetrahedra that are connected by two 
of their corners to form long chains and/or rings^ in 
agreement with experimental findings In previous 

work^ we studied the dynamical behavior of this model 
using a slightly modified set of potential parameters. We 
found that the dynamics of the atomic species decouple 
with decreasing T. While the structural relaxation of the 
phosphate network freezes in at a computer glass transi- 
tion temperature T s wl000K, fast lithium ionic diffusion 
takes place at T<T g . Thus, it is possible to equilibrate 
the lithium ionic subsystem at T <T g and to study ionic 
migration in a glassy matrix. Detailed analysis of the 
lithium dynamics showed that the dynamical behavior of 
the model resembles that of LiPC>3 glass^i 

In the present work, we identify the lithium sites in 
this model of LiPC>3 glass. For this purpose, we extend 
the temperature range of our previous study 4 ^ to lower 
T, where lithium dynamics is well described as hopping 
motion. We perform MD simulations in the NVE ensem- 
ble with p = 2.15g/cm 3 and A = 800 so that the number 
of lithium ions Nu — 160. Such moderate system size 
allows us to equilibrate the lithium ionic subsystem at 
sufficiently low T, while major finite size effects are ab- 
sent^ The density is chosen based on the experimental 
value at room temperature, p = 2.25g/cm 3 j££ together 
with the thermal expansion coefficient of LiPC>3 glass^ 
The equations of motion are integrated using the velocity 
Verlet algorithm with a time step of 2 fs, periodic bound- 
ary conditions are applied and the Coulombic forces are 
calculated via Ewald summation^ 

Our simulations start from a configuration that was 
obtained previously— by equilibrating the phosphate net- 
work slightly above T g , quenching into the glass and equi- 
librating the lithium ionic subsystem at T = 590 K. To set 
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the target temperatures of the present study, T — 390 K, 
T = 450 K and T — 520 K, we first propagate this con- 
figuration while periodically rescaling the velocities. Af- 
terwards, equilibration periods of lengths i eg = 30 — 40 ns 
without velocity scaling are applied at each T. Finally, 
we perform production runs of length t s % m = 20ns. For 
T = 520 K, all lithium ions exit their sites during the 
equilibration period t eq , indicating that the lithium ionic 
subsystem is well equilibrated. For the lower T, this con- 
dition cannot be met. However, apart from fluctuations, 
the temperature is constant after t eq and, hence, the pop- 
ulation of the lithium sites has reached its equilibrium 
distribution. Such behavior becomes reasonable based on 
our result that the number of lithium sites is very lim- 
ited, see below, so that sites with long residence times and 
low energies are occupied in any case and equilibration 
only requires the repopulation of sites with short resi- 
dence times and high energies, between which numerous 
exchange processes take place during the applied equili- 
bration periods. 

To further check the equilibration of the lithium ionic 
subsystem, we calculate the incoherent intermediate scat- 
tering functions of the lithium ions, see Fig. |3 and de- 
termine the mean time constants of the decays, as was 
done previously^ Comparison of previous and present 
results shows that all data at T <T g are well described 
by an Arrhcnius law with activation energy E a = 0.47eV. 
Fitting data for both T>T g and T<T g , we found a big- 
ger value E a = 0.62eV in our previous study^i Detailed 
inspection, however, shows that this difference can be 
traced back to a change in the temperature dependence 
at rwT g sal000K. Based on these findings, we conclude 
that the lithium ionic subsystem is equilibrated at all T 
studied here. 



III. IDENTIFICATION OF LITHIUM SITES 

Following Heuer and coworkers (2* a straightforward al- 
gorithm is applied to identify the lithium sites based on 
the MD trajectories. In detail, we first divide the simu- 
lation box into cubic cells of size (0.3 A) 3 , which is suffi- 
ciently small so as to resolve the shape of the sites. For 
each T, we then determine the number of times, m, the 
cells arc visited by a lithium ion during t s , m = 20ns, where 
configurations separated by a time increment Ai = 0.2ps 
are analyzed. The visited cells thus constitute the lithium 
sites and the pathways between them. To eliminate the 
pathways, only cells with m > m,Q are considered in a 
following cluster analysis, where cells sharing a common 
face are grouped into one site. The value of mo is de- 
termined by the criterion that the number of distinct 
lithium sites be maximum. While lithium sites that are 
occupied less frequently are not detected for high values 
of mo, different sites merge for small values of m due to 
inclusion of the connecting pathways. 

Having identified the lithium sites, the migration of 
the lithium ions can be mapped onto sequences of jumps 
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FIG. 1: Effective radius R a = (3V/47r) 1/3 of the lithium sites 
and mean number n 3 of lithium ions at a site, (a) Results 
for 178 sites obtained from a preliminary analysis at T = 
390 K. (b) Results for the final set of 166 sites determined by 
eliminating satellite and transition sites, see text for details. 
In both panels, the sites are sorted according to their effective 
radius, which is plotted in units of A, and v denotes the site 
index. 



between these sites. When an ion has left a site, it can 
either move to a new site or return to the old one. While 
we record the former event as a jump, we dismiss the 
latter and regard the residence at the particular site as 
not interrupted. In this way, we can exclude the effect 
of occasional large amplitude vibrations, which lead to 
exploration of volume outside the sites, but not to long- 
range transport. Thus, the residence time of a lithium 
ion at a site is defined as the interval between the time 
when the ion jumps into this site and the time when the 
ion exits this site in order to enter another site. Based 
on the residence times, we calculate both the mean resi- 
dence time t s at a site and the mean number of lithium 
ions at a site, or, equivalently, the occupation number n s . 
To determine r s , a sufficiently large number of jump pro- 
cesses has to be analyzed. This is possible for T = 520K, 
whereas some sites are repopulated just a few times dur- 
ing t S i m at the lower T. To reduce the effect of the limited 
time window, we multiply the residence times at the be- 
ginning and at the end of the production run by a factor 
of two when calculating t s . 

Heuer and coworkers^ observed that some small sites 
resulting from the described algorithm do not serve as 
independent sites, but they are satellites of larger sites 
or saddle-like states that are visited for a short period of 
time during the transitions between larger sites. There- 
fore, these small sites were excluded from their analysis. 
In view of these results, we determine the lithium sites 
in two steps. First, we identify the maximum number 
of sites at each T and study their nature in a prelimi- 
nary analysis. Based on the findings, we then eliminate 
satellite and transition sites and perform the final anal- 
ysis for the resulting set of sites, which turns out to be 
temperature independent. 

The preliminary analysis yields between 176 and 193 
lithium sites at the studied T. While 166 sites are present 
at all T (A sites) , there is an additional temperature de- 
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r cm < 2.0 A, which is smaller than the minimum Li-Li in- 
teratomic distance, i.e., gtiLi{f < 2.0 A) = 0, (iii) a ratio 
riij/rii < 0.3 indicates that the sites are not independently 
occupied and (iv) a value max [pij , pji] > 0.5 shows that 
jumps from one of the sites predominantly lead to the 
other. As a result, we find that all B sites identified at 
T = 390 K are associated with an A site, leading to a 
number of 166 independent lithium sites. This final set 
of lithium sites is shown in Fig. [5] We see that the sites 
are compact, mostly globular objects, which are well sep- 
arated from each other. All further analysis is based on 
these sites. However, we determined that our conclusions 
are not affected when the following analysis is performed 
on the preliminary sets of sites. Moreover, we ensured 
that the first and the second half of our production runs, 
respectively, yield nearly identical sets of sites. 



IV. RESULTS 
A. Properties of the lithium sites 

Analyzing the ionic trajectories, we find that, for all 
studied T, more than 96% of the lithium positions lie 
within the volume constituted by the 166 lithium sites, 
i.e., within a volume fraction of only 1%. These percent- 
ages show that our approach allows us a complete identi- 
fication of the lithium sites. In addition, they imply that 
the positions of the lithium sites are temperature and 
time independent, which in turn suggests that structural 
relaxation of the phosphate matrix plays no substantial 
role on the time scale of our simulation. Consistently, we 
observe that the lithium sites exhibit temperature inde- 
pendent properties, i.e., sites with relatively long (short) 
mean residence times r s , high (low) occupation numbers 
n s , etc. at a given T show the same features at the other 
T studied here. 

To investigate the dynamics of the phosphate matrix 
in more detail, we calculate the mean square displace- 
ments of the oxygen and phosphorus atoms, rg(i) and 
under the condition that site j is occupied, ny, and r ^{t), respectively Values r£ p (10ns) <0.6A 2 indicate 

that substantial structural relaxation of the phosphate 
matrix is absent at the studied T . To obtain further in- 
sights, we determine the number of phosphorus atoms 
in the first neighborshells of the oxygen atoms, i.e., we 
identify bridging oxygens (BO) and nonbridging oxygens 
(NBO). For all T, we find that more than 98% of the oxy- 
gen atoms have the same number of phosphorus neigh- 
bors at the beginning and at the end of our 20 ns-time 
window and, hence, each oxygen atom has a preferred 
phosphorus coordination. However, there are temporary 
matrix fluctuations. For example, when we analyze con- 
figurations at 1000 equidistant times during £ S j TO = 20ns, 
we observe that the fraction of oxygen atoms that never 
changes the number of phosphorus neighbors is reduced 
to 75% for T = 520K. 

Motivated by the finding that the spatial distribution 
of NBO is basically static on the time scale of our sim- 



FIG. 2: Lithium sites in LIPO3 glass. The sites were ob- 
tained by dividing the simulation box into small cubic cells 
and counting the number of times a particular cell is visited 
by a lithium ion in the course of the simulation. Frequently 
visited, face sharing cells were then grouped into one site. 
Here, the cells are shown as spheres, where cells constituting 
distinct lithium sites are given different colors. 



pendent set of sites (B sites). Specifically, when the tem- 
perature is increased, some sites merge and some new 
sites become available. Further analysis showed that 
the B sites exhibit distinctly smaller occupation numbers 
n s and effective radii R s — (3V/47T) 1 / 3 than the A sites 
so that they contribute to the drops observed for these 
quantities in Fig. ^a). The B sites are also separated 
by unusually small center-of mass distances r cm from at 
least one other site. To further study the nature of the 
lithium sites, we calculate the occupation number of site 



normalize it by the regular occupation number of site i 
Til ■ Then, a ratio n^/n, ~ 1 is observed for all A sites, in- 
dicating an independent occupation of these sites, while 
values nij/rii <C 1 are found for the B sites and, hence, 
their population is suppressed by the population of other 
sites. Finally, from inspection of the probabilities 
that a jump from site i leads to a site j, it becomes clear 
that many B sites are predominantly exited towards a 
particular other site, as is expected for satellite sites. 

These findings show that some of the lithium sites re- 
sulting from our preliminary analysis are indeed satellite 
or transition sites and, hence, their explicit consideration 
is not necessary^i To determine the "relevant" lithium 
sites, we start from the set identified at T = 390 K and 
combine any sites i and j if at least three of the following 
four criteria are fulfilled: (i) the sites merge at higher T, 
(ii) the sites are separated by a center-of mass distance 
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lgx s , lgt[ps] 

FIG. 3: Solid symbols: Cumulative distribution C(lgr s ) char- 
acterizing the mean residence times of the lithium ions at the 
sites. Open symbols: Incoherent intermediate scattering func- 
tion of the lithium ions, F a (lgt), for g = 2 7 r/r LlLl = 2.3A- 1 . 



ulation, we next study the relative positions of lithium 
sites and NBO. For this purpose, we identify the oxygen 
atoms in the hrst ncighborshells of lithium ions residing 
at a site and determine the fraction of NBO among this 
subset of oxygens. A statistical analysis shows that the 
fraction amounts to more than 96% at the studied T and, 
hence, a lithium ion at a site is predominantly surrounded 
by NBO. In other words, the lithium sites are located 
in regions with high local concentrations of NBO. For 
comparison, as expected for the meta-phosphate com- 
position, the fraction of NBO among all oxygens in the 
system amounts to 2/3. 

Figure ^b) displays the effective radii R s of the 166 
lithium sites together with their occupation numbers n s 
at T — 390 K. We see that most sites are characterized 
by radii i? s «0.5A, indicating that the sites are compact 
objects of similar size, consistent with the appearance of 
Fig. |21 Further, we observe that the vast majority of the 
sites exhibits occupation numbers n s ^ 1, i.e., they can 
accommodate only one ion at a time, so that the total 
number of available sites, N s , is approximately given by 
the number of identified sites. Hence, the number of 
sites, N s w 166, is only slightly bigger than the number 
of ions, A^Li = 160. Consistently, Heuer and coworkers^ 
reported N s w A^i for Li 2 Si03 glass and proposed to 
describe lithium dynamics in terms of mobile vacancies. 
Moreover, Dyre^ argued based on general reasoning that 
a minimum number of alkali sites is formed during the 
freezing of the network at T g due to energetic reasons. 

Next, we characterize the mean residence times r s at 
the lithium sites by the cumulative distribution C(lgr s ). 
The latter is related with the distribution of mean resi- 
dence times, G(lgT s ), by 

C(lgr,)= I 5 *' G(lgr' s ) d\gr' a 

J — oo 



so that the cumulative distribution measures the fraction 
of sites with residence times smaller than t s . In Fig. [3] 
we see that C(lg r s ) shifts to longer times with decreasing 
T, reflecting the slowing down of lithium jump dynamics. 
Furthermore, the rise of C(lgr s ) is more gradual at lower 
T and, hence, the distribution G(lgr s ) broadens upon 
cooling. In particular, the increase of C(lgr s ) extends 
over more than four orders of magnitude at T — 390 K, 
indicating that the sites feature very diverse lithium jump 
dynamics, i.e., pronounced dynamic heterogeneities exist. 
In addition, it becomes clear from Fig.[3]that about 20% 
of the lithium ions uninterruptedly reside at the same site 
during our production run for T = 390 K, and, thus, the 
mean residence time of these sites cannot be determined. 
In the following analysis, we use r s =40 ns = 2t S £ m , which 
may somewhat affect the results for T — 390 K. 

For comparison, we consider the incoherent intermedi- 
ate scattering function of the lithium ions, 

F q {t) = {cos{q[r(t +t)-r(t )}}). 

Here, r(t) is the position of a lithium ion at a time t and 
the brackets (...) denote the ensemble average. To study 
lithium dynamics on the length scale of the Li-Li inter- 
atomic distance, rnLi = 2.1 'A, we use an absolute value of 
the wave vector q = 2ir/rLiLi, as was done previously^ In 
Fig. we show the temperature dependence of F q (\gt). 
Comparing C(lgr s ) and F q (\gt) for T = 450 K, we see 
that the latter function does not decay to zero on a time 
scale of 10 ns (lgffps] =4) and, hence, a significant frac- 
tion of lithium ions still occupies the initial site, whereas 
C(lgr s [ps] = 4)k1 indicates that the residence times at 
nearly all sites were interrupted by jumps to other sites 
during this time window. In combination, these findings 
suggest that there is a substantial fraction of unsuccessful 
escape processes, i.e., the ions jump back to the initial 
site. Thus, the results of this section confirm our pre- 
vious findings^ that both dynamic heterogeneities and 
correlated back-and-forth jumps are important features 
of lithium jump dynamics in LiP03 glass. As will be 
demonstrated below, the present approach allows us to 
study the origin of these phenomena in detail. 



B. Mechanism for the lithium jumps 

To analyze the mechanism for the lithium jumps, we 
determine the mean jump time Tj and the mean jump 
length lj, where we define the former as the time an ion 
spends outside any site when it moves from one site to 
another and the latter as the center-of mass distance of 
two consecutively visited sites. At the studied T, we find 
mean jump times Tj =0.6 — 1.7 ps that are much shorter 
than the mean residence times r s , cf. Fig. |21 confirm- 
ing that jump diffusion occurs. The mean jump lengths 
lj = 3.2 — 3. 7 A are only slightly bigger than the inter- 
atomic distance rLiLi = 2.7A and, thus, the lithium ions 
predominantly jump to neighboring sites. 
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FIG. 4: Mean coordination numbers (a) 2Lio an d (b) ZLiP for 
lithium ions residing at sites with short (S T a = 1), medium 
(S T s = 2) and long (S T3 = 3) mean residence times t s . Each 
subset contains about 1/3 of all sites. For example, S T s = 1 
consist of sites with residence times t s <ti, where C(lg n) = 
1/3. 

When we consider the very limited number of lithium 
sites N s ~Nu, there are two possible mechanisms for the 
repopulation of the sites. On the one hand, one can imag- 
ine a vacancy mechanism where a site is first exited by 
ion A and immediately afterwards occupied by ion B. In- 
deed, studying examples of ionic trajectories, a vacancy- 
like mechanism was observed for a silicate glass^i On the 
other hand, it is possible that ion B enters an occupied 
site and "kicks out" ion A. Both scenarios can be distin- 
guished when we determine the probability that a lithium 
jump leads to an empty site. We find that this proba- 
bility is higher than 84% in the considered temperature 
range, indicating that a vacancy mechanism dominates 
the repopulation of the lithium sites. 

C. Origin of the dynamic heterogeneities 

The relation between the lithium jump dynamics and 
the local glass structure can be studied when we analyze 
the environments of sites with diverse mean residence 
times. Therefore, we now compare the mean coordina- 
tion numbers of lithium ions at sites with short (S TS = 1), 
medium (S TS — 2) and long (S TS = 3) residence times t s . 
In Fig. 01 we sec that, on average, lithium ions at sites 
with short r s exhibit smaller mean coordination num- 
bers ZLiO and ZLiP than those ions at sites with long r s . 
This suggests that sites featuring slow lithium jumps are 
well embedded in the phosphate-glass matrix, however, 
the effects are weak. For the mean coordination number 
Zhihi, we observe no significant correlation. Also, we find 
no systematic effects due to other quantities character- 
izing the local lithium environments. Hence, our results 
give no evidence that any particular structural property 
determines the lithium jump dynamics in L1PO3 glass. 
However, the local structure determines the local energy 



FIG. 5: Mean residence times r s at the lithium sites plotted 
as a function of (a) the site energy e s and (b) the average site 
energy of the neighboring sites, (e„). 



landscape. We find that the average energy of a lithium 
ion at a site, or, equivalently, the site energy e s is rel- 
atively low for sites where the lithium ions have oxygen 
neighbors that contain a high fraction of NBO. In the fol- 
lowing, we directly study the effects of the local energy 
landscape on lithium jump dynamics. 

To determine the site energy e s of a given site and for 
a given T, we identify the lithium ion, if any, that oc- 
cupies this site in a considered MD configuration at this 
T, calculate its potential energy using Eq. ^ and aver- 
age over many configurations. Figure [BT a) shows t s as a 
function of the site energies e s . As may be expected, the 
mean residence time decreases with increasing e s . How- 
ever, decays by factors of between five and ten are small 
compared to the width of G(\gr s ), which spans at least 
three orders of magnitude, see Fig. [3J For example, if 
a random trap model were an appropriate description, 
i.e., if the residence time were exclusively determined by 
e s , a decline of r s by a factor of ~ 2000 would be ex- 
pected based on Ar s = exp[Ae s /(ksT)], Ae s = 0.3 eV 
and T = 450 K. Thus, the jump rate at a site depends 
not only on e s , but there are additional factors. For ex- 
ample, the energy barriers separating the sites may show 
a distribution of barrier heights, as expected for the en- 
ergy landscape of a disordered material. 

Another important factor can be inferred from an in- 
spection of Fig. GJb) , where the mean residence time at 
a site is plotted as a function of the average site energy, 
of all neighboring sites. Here, we define neighboring 
sites based on the criterion that their center-of mass dis- 
tance r cm be smaller than the distance corresponding to 
the first minimum of the Li-Li pair correlation function. 
Interestingly, we observe that t s strongly decreases with 
increasing (e„), indicating that, on average, sites sur- 
rounded by high energetic neighbors feature fast lithium 
dynamics. On first glance, this result is counterintuitive 
since one would expect that jumps to sites with compar- 
atively high site energies are characterized by small rates 
and, hence, sites with high values of (e„) should exhibit 
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FIG. 6: (a) Mean number n s of lithium ions at a site as a 
function of the site energy e a . (b) Mean residence time r a at 
a lithium site as a function of n s . 

long r a . 

However, a decrease of t s with can be understood 
when we consider our result that a vacancy mechanism 
dominates the redistribution of a small fraction of empty 
sites. In such scenario, the ionic jump rates are first of all 
limited by the availability of a vacant neighboring site. 
In Fig.[Ua), we see that the occupation numbers n s are 
smaller for sites with high site energies s s , as may be 
expected. Thus, when the neighbors of a site A exhibit 
high site energies, these neighbors are often unoccupied 
so that ions at site A frequently have the chance for a 
jump. Consequently, in our case of a vacancy mecha- 
nism, it is plausible that sites characterized by high val- 
ues of (e n ) show relatively short r s , as is observed in 
Fig- Elk)- To further check our argumentation, we plot 
the mean residence time at a site as a function of the av- 
erage occupation number, (n n ), of all neighboring sites in 
Fig-EHb). We see that t s strongly grows with (n n ), sup- 
porting our conclusion that lithium diffusion in LiPC>3 
glass is governed by the competition of the ions for a 
small fraction of vacant sites. These results clearly show 
that single-particle approaches do not allow a complete 
understanding of ion dynamics in glasses. 

D. Correlated back-and-forth jumps 

The results of our previous worki^ and the discus- 
sion of Fig. pimply that correlated back-and-forth jumps 
are another feature of lithium dynamics in the studied 
model. To investigate this effect in more detail, we cal- 
culate the backjump probability P of finding a lithium 
ion at the same site after exactly two jumps. Analyz- 
ing all jumps during our production runs, we find that 
the backjump probability increases from P b = 0.77 at 
T = 520 K to P b = 0.91 at T = 390 K. For comparison, 
a value P b rs 0.25 would result for a random walk since, 
on average, each lithium site has about four neighboring 
sites. Thus, the observed backjump probabilities clearly 



FIG. 7: Backjump probability p b as a function of the loga- 
rithm of the mean residence time, lgr s , for various temper- 
atures. Dashed line: Cumulative distribution of the mean 
residence time, C(lgr s ), for T = 520K. 

indicate that correlated back-and-forth jumps are an im- 
portant phenomenon in the studied temperature range. 

To obtain further insights, we next determine the site- 
dependent backjump probability p b , i.e., for each site, 
we calculate the probability that a jump from this site to 
another site is followed by a direct backjump. Motivated 
by results of Heuer and coworkers^ we then study the 
relation between the backjump probability p b and the 
mean residence time t s . From the results in Fig. [7J it is 
evident that the backjump probability is higher at lower 
T. Furthermore, we see that pjj strongly decreases with 
increasing lg t s . Specifically, comparison with the cumu- 
lative distribution C(lgr s ) shows that about 90% of the 
jumps starting from sites with short lgr s are followed by 
a direct backjump, whereas backjumps occur with nearly 
statistical probability, p b « 0.25, after the escape from 
sites with long r s . Reinspecting also F q (\gt) in Fig. |3J 
these findings show that lithium hopping on the time 
scale of the lithium relaxation resembles a random walk, 
whereas ions that exit their sites at much earlier times 
very often return to these sites, i.e., such jumps are un- 
successful. A similar behavior was observed for lithium 
silicate glasses by Heuer and coworkers & 

Finally, we analyze the relation between the backjump 
probability and the local energy landscape. In Fig.|Sfa), 
we see that, on average, sites with high site energies e s 
show somewhat higher values of p b . Figure |8fb) displays 
p b as a function of the average energy, (e n ), of the neigh- 
boring sites. Obviously, the backjump probability hardly 
depends on (e n ). When we assume single-particle motion 
in an energy landscape, these results are difficult to un- 
derstand. In such approach, jumps to sites with lower 
site energies are favored and, hence, the backjump prob- 
ability should be high for sites characterized by small 
values of e s and high values of (e n ), at variance with our 
observations. Hence, our findings for the back-and-forth 
dynamics again suggest that multiparticle interactions 
are important for the lithium jumps in LiP03 glass. We 
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FIG. 8: Backjump probability p h a as a function of (a) the site 
energy e s and (b) the average site energy of the neighboring 
sites, (e n ). 



add that we find no evidence for a systematic dependence 
of Pg on properties of the local glass structure such as the 
mean coordination numbers ZLiO: z LiP and ZLiLi- 



V. DISCUSSION AND SUMMARY 

We performed MD simulations to investigate lithium 
dynamics in L1PO3 glass at temperatures where the phos- 
phate matrix is basically rigid on the time scale of our 
simulation. We observed that lithium ionic migration in 
this glassy matrix results from jumps between well de- 
fined sites. Making use of a straightforward algorithm, 51 
the lithium sites were completely identified from the ionic 
trajectories. The results show that the sites are typically 
compact objects of similar size and globular shape. More- 
over, the number of lithium sites hardly outnumbers the 
number of lithium ions, as was observed for lithium sili- 
cate glasses^ Characterizing the structural and dynam- 
ical features of the sites, we found that these properties 
are unmodified on the time scale of the lithium ionic re- 
laxation. This allowed us to systematically study the 
relationship between structure and dynamics. 

Analyzing the phosphate dynamics in detail, we found 
that the glassy matrix, apart from temporary fluctu- 
ations, is static on the 20ns-time scale of our simu- 
lation. In particular, there is a basically time inde- 
pendent spatial distribution of NBO. While one may 
speculate whether the temporary matrix fluctuations as- 
sist the lithium jumps, as was observed for other ionic 
glassesjSSiSi the present results are at variance with a sub- 
stantial structural relaxation of the glassy matrix, e.g., 
a redistribution of NBO, as was proposed in previous 
workpiSiiS Instead, the properties of the lithium sites are 
determined by the configuration frozen in at the glass 
transition. In particular, we observed that the lithium 
sites are located in regions with high local concentra- 



tions of NBO so that most lithium ions at a site are 
exclusively surrounded by this oxygen species. Though 
we do not expect that our findings are unique to LiP03 
glass, we cannot exclude that the scenario changes for 
small alkali concentrations, where NBO may be formed 
along the alkali trajectories 

Mapping the lithium trajectories onto sequences of 
jumps between the identified sites, we analyzed the mech- 
anism for the ionic diffusion. The mean jump times and 
the mean jump lengths indicated that lithium migration 
results from hopping between neighboring sites. A sta- 
tistical analysis of the repopulation of the lithium sites 
showed that the sites are usually first vacated by one ion 
and shortly afterwards entered by another ion, implying 
that lithium hopping is governed by a vacancy mecha- 
nism. Such mechanism is reasonable in a situation where 
the number of alkali sites is only slightly bigger than the 
the number of alkali ions, as is the case both in lithium 
silicate glasses^ and in our model of a lithium phosphate 
glass. For silicate glasses, a vacancy-like mechanism was 
also observed when inspecting examples of ionic trajec- 
tories^ 

To characterize the average dynamical behavior of a 
lithium ion at a particular site, we determined the mean 
residence times r s at the sites. We found that the sites 
feature a broad distribution G(lgr s ), indicating the ex- 
istence of pronounced dynamic heterogeneities. G(lgr s ) 
broadens upon cooling so that the distribution extends 
over more than four orders of magnitude at the lowest 
T studied. Consistently, we observed in our previous 
MD study of LiP03 glass^ that dynamic heterogeneities 
increasingly contribute to the nonexponential lithium re- 
laxation upon cooling. All these simulation results are 
in agreement with findings from multidimensional NMR 
experiments on silver phosphate glassesjSL^i which show 
that a very broad distribution of jump rates governs sil- 
ver dynamics at low T, where the silver relaxation occurs 
on the time scale of ms — s. 

The present approach provides direct access to the ori- 
gin of the dynamic heterogeneities. We observed that, 
on average, lithium ions at sites with long r s show some- 
what higher coordination numbers zuo and ^LiP, sug- 
gesting that these sites are located in a matrix-rich en- 
vironment. However, the direct effects of the local glass 
structure on the lithium jump dynamics are weak. On 
the other hand, the local structure determines the local 
energy landscape, which in turn strongly affects the res- 
idence times. Specifically, as may be expected, we found 
that r s decreases with the site energy e s . More interest- 
ingly, t s also declines when the average site energy (e n ) 
of the neighboring sites increases. Such behavior is dif- 
ficult to capture within single-particle approaches on ion 
transport in glasses, including percolation approaches. 
However, it can be understood when we consider our re- 
sults that there is a small fraction of empty sites that are 
redistributed in a vacancy mechanism. In such scenario, 
the jump rates are limited by the restricted access to a 
vacant site. Due to a correlation between the energies e s 
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and the occupation numbers n s of the sites, the proba- 
bility of finding a vacant neighboring site is high when 
the latter exhibit high site energies and, consequently, t s 
is expected to decrease with (e n ), in agreement with our 
observation. We conclude that the pronounced dynamic 
heterogeneities in LiPC>3 glass result from both a disor- 
dered energy landscape and a distribution of probabilities 
of having access to a vacant neighboring site. 

We also calculated the site-dependent backjump prob- 
ability p s b that a jump from a particular site to another 
site is followed by a direct backjump. Such analysis al- 
lowed us to quantify the importance of correlated back- 
and-forth jumps. We found that p% strongly depends on 
the mean residence time r s , i.e., on the dynamical state of 
an ion. High values of p% for sites with short r s indicate 
that correlated back-and-forth jumps are a relevant fea- 
ture of lithium dynamics at the studied T. Specifically, 
backjump probabilities p s b ~0.9 that increase upon cool- 
ing show that only ~ 10% of the jumps starting from sites 
with short r s are successful and, hence, ion transport is 
substantially slowed down due to back-and-forth jumps, 
in particular at low T. In contrast, sites characterized 
by long t s show backjump probabilities that are close to 
the value expected for a random jump. These results for 
LiPC>3 glass are in agreement with findings for lithium 
silicate glasses by Heuer and coworkers j** implying that 
there is a similar mechanism for lithium jump dynamics 
in different types of glasses. However, correlated back- 
and-forth jumps were not observed in a computational 
study of ion dynamics at comparatively high T4£ This 
apparent discrepancy may be reconciled by our result 
that the importance of this phenomenon decreases with 
increasing T. 

Our results may also help to unravel two puzzling ef- 
fects related to ion dynamics. As aforementioned, appli- 
cations of different experimental techniques lead to dif- 
ferent conclusions about the relevance of correlated back- 
and-forth jumps. These discrepancies may be explained 
by the strong dependence of the backjump probability 
on the dynamical state of an ion. Similar to Heuer and 
coworkers^ we argue that the observed correlated back- 
and-forth jumps of the fast ions dominate the dispersive 
regime of the electric conductivity ct(V)i 6 i 24 i 25 i 26 i 27 An 
additional contribution to the dispersion may result from 
fast forward-backward motions on length scales shorter 
than the Li-Li interatomic distance, e.g., jumps between 
sites and their satellites, which are not considered in the 
present analysis. In contrast, studying the dynamical be- 
havior of subensembles of slow ions in multidimensional 
NMR experiments no evidence for correlated back- 



and-forth jumps was found. Consistently, we observed 
that lithium dynamics on the time scale of the lithium 
relaxation resembles a random walk. 

The second puzzling effect concerns the approximate 
time-temperature superposition observed for the incoher- 
ent intermediate scattering function F q (t) in simulation 
studies42ii4 This behavior is surprising since both dy- 
namic heterogeneities^ii and back-and-forth jumps be- 
come more important with decreasing T, see Figs. Eland 
and, hence, one might be tempted to conclude that 
the nonexponentiality of F q (t) should increase upon cool- 
ing. However, such simple argumentation is only possible 
when the backjump probability is independent of the dy- 
namical state. Here, the back-and-forth jumps of the 
fast ions shift the time to when F q (t) starts to decay to 
longer times, while the time t\ when the decay of F q (t) is 
complete is not affected due to back-and-forth dynamics 
because no enhanced back-jump probability exists for the 
slow ions. This means that, in our case, back-and-forth 
jumps reduce the time interval t\— to, i.e., the nonexpo- 
nentiality of F q (t). Hence, when T is lowered, there are 
two competing effects. While the increasing heterogene- 
ity of dynamics leads to a stronger nonexponentiality, the 
growing backjump probability of the fast ions delays the 
onset of the decay, to, further and further and, hence, 
reduces the nonexponentiality. Thus, both effects may 
compensate each other, resulting in a nearly tempera- 
ture independent nonexponentiality of F q (t). 

In summary, we found that, in a model of LiPC>3 
glass, the disordered energy landscape provided by the 
phosphate-glass matrix plays a major role for the lithium 
jump dynamics. For a complete understanding of the 
jump-diffusion mechanism, single-particle approaches are 
not sufficient, but rather it is important to consider the 
competition of the lithium ions for a small fraction of 
vacant sites at every instant. As a consequence, the dy- 
namical behavior is highly complex. In particular, pro- 
nounced dynamic heterogeneities exist and the backjump 
probability depends on the dynamical state of an ion. 
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